# delimit ; 
clear; 
set mem 300m;	
set more 1 ;  
drop _all;
program drop _all;
cd "C:\Dropbox\ERIC_couples\NewPC\Exported Sim data\";

cap log close;
log using deathwealth_rm_Simcouples5.log, replace; 

use "event_HRSTrim_230117.dta";
*use "event_HRSTrim_onebeq.dta";

drop _merge;
*First restructure Data around Events;
*To do this, it is essential that we calculate the correct wave of death-wave;
*Because households actually die between waves, use 2nd wave;

gen manwave=wave if mandied==1;
gen womanwave=wave if womandied==1;

*use the fact that max of missings and nonmissing returns the non missing value;
egen mandiedwave=max(manwave), by(HHID);
egen womandiedwave=max(womanwave), by(HHID);

egen firstwave=min(wave), by(HHID);

*generate evercouple id;

drop couple;
gen couple=(hhstatus==3);

egen evercouple=max(couple), by(HHID);

*The first and last death (only really relevant for couples);
gen firstdiedwave=min(womandiedwave,mandiedwave) if evercouple==1;
gen lastdiedwave=max(womandiedwave,mandiedwave);

*Households who are "single", but first appear at a spousal death;
replace lastdiedwave=. if lastdiedwave==firstwave & evercouple==0;

* drop manwave womanwave;

*Now make sure that lastdied and firstdied only apply appropriately;
*replace if they are picking up the same event i.e. lastdied hasn't happened yet;
gen indicator=(hhstatus==0);

egen observe_alldead=max(indicator), by(HHID);
replace lastdiedwave=. if evercouple==1 & observe_alldead==0; 

drop observe_alldead indicator;

*restructure around mandied and womandied;
*These variables all contain the calender time dimension re
gen event_mandied=wave-mandiedwave;
gen event_womandied=wave-womandiedwave;
gen event_lastdied=wave-lastdiedwave;
gen event_firstdied=wave-firstdiedwave;

rename OOP hhmedcost_b;
rename OOP_don hhmedcost_b_don;
rename Medicaid hhmedicaid_pay;
rename Medicaid_don hhmedicaid_pay_don;

* correct concept is household out of pocket + Medicaid;
rename hhmedcost_b medcost_b;
rename hhmedicaid_pay medicaid_pay;

rename hhmedcost_b_don medcost_b_don;
rename hhmedicaid_pay_don medicaid_pay_don;

*select the dependent variable;
gen medcost=medcost_b+medicaid_pay;
gen medcost_don=medcost_b_don+medicaid_pay_don;

gen years_last=2*event_lastdied;
gen years_first=2*event_firstdied;

la var years_last "Years";
la var years_first "Years";

*Adjust for the way in which we record death assets (pre beq choice) in simulations;
replace assets=assets-beq if wave!=lastdiedwave;
*updated 230117;


rename pi_perc PI;
*keep HHID wave asset* homee* medcos* medicaid* manwave womanwave event_* *diedwave years_last years_first evercouple PI;
keep HHID* wave asset* medcos* medicaid* manwave womanwave event_* *diedwave years_last years_first evercouple PI amothers* age;



*expand 2, gen(expander) ;

*la def expandlab 1 "No Death" 0 "Death in Household";

*lab val expander expandlab;
/*
replace assets=assets_don if expander==1;
replace medcost=medcost_don if expander==1;
replace medcost_b=medcost_b_don if expander==1;
replace medicaid_pay=medicaid_pay_don if expander==1;
replace homeeq=homeeq_don if expander==1;
*/
*replace amothers=amothers_don if expander==1;

*replace amothers=amothers/1000;
*la var amothers "Amount Not Transfered to Spouse at Death (1000s 2014 Dollars)";


gen w3wealth = assets if wave ==3;
egen w3iwealth = max(w3wealth), by(HHID);
drop w3wealth;

gen deathage1 = age if years_first==0 & evercouple==1;
drop age;
egen deathage = max(deathage1), by(HHID);

gen age=deathage;

egen median_assets_couples=median(w3iwealth) if evercouple==1;
gen below_median_assets_couples=1 if w3iwealth<median_assets_couples & evercouple==1;
gen above_median_assets_couples=1 if w3iwealth>median_assets_couples & evercouple==1;

egen median_assets_singles=median(w3iwealth) if evercouple==0;
gen below_median_assets_singles=1 if w3iwealth<median_assets_singles & evercouple==0;
gen above_median_assets_singles=1 if w3iwealth>median_assets_singles & evercouple==0;



*replace homeeq=homeeq/1000;
*replace homeeq_don=homeeq_don/1000;

*gen finwealth=assets-homeeq;
*gen finwealth_don=assets_don-homeeq_don;

*****************************************************************;
* estimate wealth decline\medical expense rise right before death;
*****************************************************************;
* robustness checks                                              ;
* top coding for assets                                          ;
* what to do with households where mstat=3 in wave 1             ;
* keep wave 1, versus drop it                                    ;
*****************************************************************;

* we have some problems with missings for wave 6.  get rid of variables where the individual died in wave 6, and wave 6 data are missing;

sort HHID wave;
*keep if HHID < 200470;

*replace assets*=. if years_first <-6;
*replace assets*=. if years_first > 4;

*drop if years_first <-6;
*drop if years_first > 4;


* generate lagged death (periods before death), single -->dead;
/*
gen singlediedtmo=0; * tmo = "t-1"; 
replace singlediedtmo=1 if years_last[_n+2]==0 & evercouple==0 & HHID==HHID[_n+2];

gen singlediedtmt=0; * tmt = "t-2"; 
replace singlediedtmt=1 if singlediedtmo[_n+2]==1 & evercouple==0 & HHID==HHID[_n+2] ;

gen singlediedtmth=0; * tmth = "t-3"; 
replace singlediedtmth=1 if singlediedtmt[_n+2]==1 & evercouple==0 & HHID==HHID[_n+2];

gen singlediedbefore=0; *waves before "t-3";
replace singlediedbefore=1 if singlediedtmth[_n+2]==1 & evercouple==0 & HHID==HHID[_n+2];
replace singlediedbefore=1 if singlediedbefore[_n+2]==1 & evercouple==0 & HHID==HHID[_n+2];
*/

gen singlediedtmo=0; * tmo = "t-1"; 
replace singlediedtmo=1 if years_last==-2 & evercouple==0;

gen singlediedtmt=0; * tmt = "t-2"; 
replace singlediedtmt=1 if years_last==-4 & evercouple==0;

gen singlediedtmth=0; * tmth = "t-3"; 
replace singlediedtmth=1 if years_last==-6 & evercouple==0;

gen singlediedbefore=0; *waves before "t-3";
replace singlediedbefore=1 if years_last<-6 & evercouple==0;

gen singledied=0;
replace singledied=1 if years_last==0 & evercouple==0;

gen singleafter=0; *waves after "t";
replace singleafter=1 if years_last>0 & evercouple==0;


/*
* generate lagged death (periods before death), couple -->single;
gen onediedtmo=0; * tmo = "t-1"; 
replace onediedtmo=1 if years_first[_n+2]==0 & evercouple==1 & HHID==HHID[_n+2];

gen onediedtmt=0; * tmt = "t-2"; 
replace onediedtmt=1 if onediedtmo[_n+2]==1 & evercouple==1 & HHID==HHID[_n+2] ;

gen onediedtmth=0; * tmth = "t-3"; 
replace onediedtmth=1 if onediedtmt[_n+2]==1 & evercouple==1 & HHID==HHID[_n+2];

gen onediedbefore=0; *waves before "t-3";
replace onediedbefore=1 if onediedtmth[_n+2]==1 & evercouple==1 & HHID==HHID[_n+2];
replace onediedbefore=1 if onediedbefore[_n+2]==1 & evercouple==1 & HHID==HHID[_n+2];
*/

gen onediedtmo=0; * tmo = "t-1"; 
replace onediedtmo=1 if years_first==-2 & evercouple==1;

gen onediedtmt=0; * tmt = "t-2"; 
replace onediedtmt=1 if years_first==-4 & evercouple==1;

gen onediedtmth=0; * tmth = "t-3"; 
replace onediedtmth=1 if years_first==-6 & evercouple==1;

gen onediedbefore=0; *waves before "t-3";
replace onediedbefore=1 if years_first<-6 & evercouple==1;

gen onedied=0;
replace onedied=1 if evercouple==1 & years_first==0;

* generate periods after death, couple -->single;
/*
sort HHID wave;
gen onediedtpo=0; * tpo = "t+1"; 
replace onediedtpo=1 if years_first[_n-2]==0 & evercouple==1 & HHID==HHID[_n-2];

gen onediedtpt=0; * tpt = "t+2"; 
replace onediedtpt=1 if onediedtpo[_n-2]==1 & evercouple==1 & HHID==HHID[_n-2] ;

gen onediedafter=0; *waves after "t+2";
replace onediedafter=1 if onediedtpt[_n-2]==1 & evercouple==1 & HHID==HHID[_n-2];
replace onediedafter=1 if onediedafter[_n-2]==1 & evercouple==1 & HHID==HHID[_n-2];
*/

sort HHID wave;
gen onediedtpo=0; * tpo = "t+1"; 
replace onediedtpo=1 if years_first==2 & evercouple==1;

gen onediedtpt=0; * tpt = "t+2"; 
replace onediedtpt=1 if years_first==4 & evercouple==1;

gen onediedafter=0; *waves after "t+2";
replace onediedafter=1 if years_first>4 & evercouple==1;

sum singledied* onedied*;

sort wave;
by wave: sum singledied* onedied*;

sort HHID wave;



*divides by 2 annualize medical expenses -- this is already done in the data from the HRS;

replace medcost = medcost/2;
replace medcost_don = medcost_don/2;
la var medcost "Medical Expenses (1000s of 2014 Dollars)";
replace medcost_b = medcost_b/2;
replace medicaid_pay = medicaid_pay/2;
replace medcost_b_don = medcost_b_don/2;
replace medicaid_pay_don = medicaid_pay_don/2;

****************************************************************************************************;
*                                                     asset regressions for wave 2 singles who died ;
****************************************************************************************************;


*gen singledied=1 if evercouple==0 & hhstatus==0 & (hhstatus[_n-2]==1|hhstatus[_n-2]==2) & HHID==HHID[_n-2];


*gen onedied = 1 if evercouple==1 & (hhstatus==1 | hhstatus==2) & hhstatus[_n-2]==3 & HHID==HHID[_n-2];


replace assets=. if singleafter==1;
replace assets_don=. if singleafter==1;
replace assets=. if onediedafter==1;
replace assets_don=. if onediedafter==1;

replace assets=. if onediedbefore==1;
replace assets_don=. if onediedbefore==1;

replace assets=. if singlediedbefore==1;
replace assets_don=. if singlediedbefore==1;


replace medcost=. if onediedbefore==1;
replace medcost_don=. if onediedbefore==1;


replace medcost=. if singleafter==1;
replace medcost_don=. if singleafter==1;
replace medcost=. if onediedafter==1;
replace medcost_don=. if onediedafter==1;

replace medcost=. if singlediedbefore==1;
replace medcost_don=. if singlediedbefore==1;

replace medcost_b=. if singleafter==1;
replace medcost_b_don=. if singleafter==1;

replace medcost_b=. if onediedafter==1;
replace medcost_b_don=. if onediedafter==1;

replace medcost_b=. if onediedbefore==1;
replace medcost_b_don=. if onediedbefore==1;

replace medcost_b=. if singlediedbefore==1;
replace medcost_b_don=. if singlediedbefore==1;

replace medicaid_pay=. if singleafter==1;
replace medicaid_pay_don=. if singleafter==1;

replace medicaid_pay=. if onediedafter==1;
replace medicaid_pay_don=. if onediedafter==1;

replace medicaid_pay=. if onediedbefore==1;
replace medicaid_pay_don=. if onediedbefore==1;

replace medicaid_pay=. if singlediedbefore==1;
replace medicaid_pay_don=. if singlediedbefore==1;



replace amothers=. if singleafter==1;
replace amothers_don=. if singleafter==1;

replace amothers=. if onediedafter==1;
replace amothers_don=. if onediedafter==1;

replace amothers=. if onediedbefore==1;
replace amothers_don=. if onediedbefore==1;

replace amothers=. if singlediedbefore==1;
replace amothers_don=. if singlediedbefore==1;




/*
replace assets_don=. if expander==0;
replace assets=. if expander==1;
replace medcost_don=. if expander==0;
replace medcost=. if expander==1;
*/
/*
replace assets_don=. if assets_don!=assets;
replace assets=. if assets==assets_don;
replace medcost_don=. if medcost_don!=medcost;
replace medcost=. if medcost==medcost_don;
*/

set scheme s1mono;

set scheme s1mono;

* fixed effects;
local iter = 1;
foreach varname in "assets" "assets_don" "medcost" "medcost_don" "medcost_b" "medcost_b_don" "medicaid_pay" "medicaid_pay_don" {;

sort wave;
by wave: sum `varname';

sort HHID;

*replace prob6=1 if hhstatus==0 & wave==6 & `varname'==.;
* also, make the comparrison variable missing if  
*replace prob6=1 if hhstatus==0 & wave==6 & `iter'==2 & assets==.;
*replace prob6=1 if hhstatus==0 & wave==6 & `iter'==4 & medcost==.;

*egen problem`iter'=sum(prob6), by(HHID);

/*
if use_wave_6==1{;
replace `varname'=. if problem`iter'~=0 & alive5==1;
};
if use_wave_6==0{;
replace `varname'=. if  alive5==1;
};

replace `varname'=. if singlediedtmf==1;
*/
/*
xtreg `varname' singlediedtmo singlediedtmt singledied, fe i(HHID) ;
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a`iter'=.;
replace a`iter'=_b[_cons] if wave==3;
replace a`iter'=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a`iter'=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a`iter'=_b[_cons]+_b[singledied] if wave==6;
*/
*gen b`iter'=.;
*sum `varname' if onediedtmth==1 & evercouple==1, det;
*replace b`iter'=r(mean);
xtreg `varname' onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a`iter'=.;
replace a`iter'=_b[_cons] if wave==3;
replace a`iter'=_b[_cons]+_b[onediedtmt] if wave==4;
replace a`iter'=_b[_cons]+_b[onediedtmo] if wave==5;
replace a`iter'=_b[_cons]+_b[onedied] if wave==6;
replace a`iter'=_b[_cons]+_b[onediedtpo] if wave==7;
replace a`iter'=_b[_cons]+_b[onediedtpt] if wave==8;


gen va`iter'=.;
/*replace va`iter'=_se[_cons] if wave==3;
replace va`iter'=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/

replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

****************************************************************************************************;
*                                           graphs                                                  ;
****************************************************************************************************;

if `iter' == 1 {;
gen time=wave-6;
gen kflag = 0;
gen time_before_death=time;
*replace time_before_death=. if time==-3|time==-4;
};
replace kflag = 0;
*replace time=wave-5 if wave==1|wave==2;
*replace time_before_death=time;
*replace time_before_death=. if time==-3|time==-4;
sort time HHID;
replace kflag = 1 if time==time[_n-1];


label variable time_before_death "Years from Death";


if "`varname'" == "assets_don" {;

label variable a1 "Death in Household";
label variable a2 "No Death";



};



display "`iter'";


local iter = `iter' + 1;
};


*gen above_assets=.;
*sum assets if onediedtmth==1 & above_median_assets_couples==1, det;
*replace above_assets=r(mean);
xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if above_median_assets_couples==1, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a13=.;
replace a13=_b[_cons] if wave==3;
replace a13=_b[_cons]+_b[onediedtmt] if wave==4;
replace a13=_b[_cons]+_b[onediedtmo] if wave==5;
replace a13=_b[_cons]+_b[onedied] if wave==6;
replace a13=_b[_cons]+_b[onediedtpo] if wave==7;
replace a13=_b[_cons]+_b[onediedtpt] if wave==8;


gen va13=.;
/*replace va13=_se[_cons] if wave==3;
replace va13=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va13=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va13=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va13=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va13=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;

local iter=13;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;


sort wave;
by wave: sum assets if above_median_assets_couples==1;

*gen above_assets_don=.;
*sum assets_don if onediedtmth==1 & above_median_assets_couples==1, det;
*replace above_assets_don=r(mean);
xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if above_median_assets_couples==1, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a14=.;
replace a14=_b[_cons] if wave==3;
replace a14=_b[_cons]+_b[onediedtmt] if wave==4;
replace a14=_b[_cons]+_b[onediedtmo] if wave==5;
replace a14=_b[_cons]+_b[onedied] if wave==6;
replace a14=_b[_cons]+_b[onediedtpo] if wave==7;
replace a14=_b[_cons]+_b[onediedtpt] if wave==8;


gen va14=.;
/*replace va14=_se[_cons] if wave==3;
replace va14=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va14=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va14=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va14=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va14=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=14;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if above_median_assets_couples==1;

*gen below_assets=.;
*sum assets if onediedtmth==1 & below_median_assets_couples==1, det;
*replace below_assets=r(mean);
xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if below_median_assets_couples==1, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a15=.;
replace a15=_b[_cons] if wave==3;
replace a15=_b[_cons]+_b[onediedtmt] if wave==4;
replace a15=_b[_cons]+_b[onediedtmo] if wave==5;
replace a15=_b[_cons]+_b[onedied] if wave==6;
replace a15=_b[_cons]+_b[onediedtpo] if wave==7;
replace a15=_b[_cons]+_b[onediedtpt] if wave==8;

gen va15=.;
/*replace va15=_se[_cons] if wave==3;
replace va15=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va15=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va15=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va15=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va15=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=15;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

by wave: sum assets if below_median_assets_couples==1;

*gen below_assets_don=.;
*sum assets_don if onediedtmth==1 & below_median_assets_couples==1, det;
*replace below_assets_don=r(mean);
xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if below_median_assets_couples==1, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a16=.;
replace a16=_b[_cons] if wave==3;
replace a16=_b[_cons]+_b[onediedtmt] if wave==4;
replace a16=_b[_cons]+_b[onediedtmo] if wave==5;
replace a16=_b[_cons]+_b[onedied] if wave==6;
replace a16=_b[_cons]+_b[onediedtpo] if wave==7;
replace a16=_b[_cons]+_b[onediedtpt] if wave==8;

gen va16=.;
/*replace va16=_se[_cons] if wave==3;
replace va16=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va16=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va16=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va16=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va16=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=16;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if below_median_assets_couples==1;



*gen above_PI_assets=.;
*sum assets if onediedtmth==1 & evercouple==1 & PI>=0.5, det;
*replace above_PI_assets=r(mean);
xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI>=0.5, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a17=.;
replace a17=_b[_cons] if wave==3;
replace a17=_b[_cons]+_b[onediedtmt] if wave==4;
replace a17=_b[_cons]+_b[onediedtmo] if wave==5;
replace a17=_b[_cons]+_b[onedied] if wave==6;
replace a17=_b[_cons]+_b[onediedtpo] if wave==7;
replace a17=_b[_cons]+_b[onediedtpt] if wave==8;

gen va17=.;
/*replace va17=_se[_cons] if wave==3;
replace va17=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va17=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va17=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va17=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va17=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/


local iter=17;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets if evercouple==1 & PI>=0.5;

*gen above_PI_assets_don=.;
*sum assets_don if onediedtmth==1 & evercouple==1 & PI>=0.5, det;
*replace above_PI_assets_don=r(mean);
xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI>=0.5, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a18=.;
replace a18=_b[_cons] if wave==3;
replace a18=_b[_cons]+_b[onediedtmt] if wave==4;
replace a18=_b[_cons]+_b[onediedtmo] if wave==5;
replace a18=_b[_cons]+_b[onedied] if wave==6;
replace a18=_b[_cons]+_b[onediedtpo] if wave==7;
replace a18=_b[_cons]+_b[onediedtpt] if wave==8;

gen va18=.;
/*replace va18=_se[_cons] if wave==3;
replace va18=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va18=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va18=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va18=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va18=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=18;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if evercouple==1 & PI>=0.5;

*gen below_PI_assets=.;
*sum assets if onediedtmth==1 & evercouple==1 & PI<0.5, det;
*replace below_PI_assets=r(mean);
xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.5, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a19=.;
replace a19=_b[_cons] if wave==3;
replace a19=_b[_cons]+_b[onediedtmt] if wave==4;
replace a19=_b[_cons]+_b[onediedtmo] if wave==5;
replace a19=_b[_cons]+_b[onedied] if wave==6;
replace a19=_b[_cons]+_b[onediedtpo] if wave==7;
replace a19=_b[_cons]+_b[onediedtpt] if wave==8;


gen va19=.;
/*replace va19=_se[_cons] if wave==3;
replace va19=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va19=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va19=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va19=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va19=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=19;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets if evercouple==1 & PI<0.5;

*gen below_PI_assets_don=.;
*sum assets_don if onediedtmth==1 & evercouple==1 & PI<0.5, det;
*replace below_PI_assets_don=r(mean);
xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.5, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a20=.;
replace a20=_b[_cons] if wave==3;
replace a20=_b[_cons]+_b[onediedtmt] if wave==4;
replace a20=_b[_cons]+_b[onediedtmo] if wave==5;
replace a20=_b[_cons]+_b[onedied] if wave==6;
replace a20=_b[_cons]+_b[onediedtpo] if wave==7;
replace a20=_b[_cons]+_b[onediedtpt] if wave==8;


gen va20=.;
/*replace va20=_se[_cons] if wave==3;
replace va20=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va20=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va20=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va20=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va20=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=20;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if evercouple==1 & PI<0.5;

*Do this by PI tercile split;

xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.34, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a201=.;
replace a201=_b[_cons] if wave==3;
replace a201=_b[_cons]+_b[onediedtmt] if wave==4;
replace a201=_b[_cons]+_b[onediedtmo] if wave==5;
replace a201=_b[_cons]+_b[onedied] if wave==6;
replace a201=_b[_cons]+_b[onediedtpo] if wave==7;
replace a201=_b[_cons]+_b[onediedtpt] if wave==8;


gen va201=.;

local iter=201;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets if evercouple==1 & PI<0.34;

xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.34, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a202=.;
replace a202=_b[_cons] if wave==3;
replace a202=_b[_cons]+_b[onediedtmt] if wave==4;
replace a202=_b[_cons]+_b[onediedtmo] if wave==5;
replace a202=_b[_cons]+_b[onedied] if wave==6;
replace a202=_b[_cons]+_b[onediedtpo] if wave==7;
replace a202=_b[_cons]+_b[onediedtpt] if wave==8;


gen va202=.;

local iter=202;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if evercouple==1 & PI<0.34;

xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.67 & PI>=0.34, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a203=.;
replace a203=_b[_cons] if wave==3;
replace a203=_b[_cons]+_b[onediedtmt] if wave==4;
replace a203=_b[_cons]+_b[onediedtmo] if wave==5;
replace a203=_b[_cons]+_b[onedied] if wave==6;
replace a203=_b[_cons]+_b[onediedtpo] if wave==7;
replace a203=_b[_cons]+_b[onediedtpt] if wave==8;


gen va203=.;

local iter=203;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets if evercouple==1 & PI<0.67 & PI>=0.34;

xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.67 & PI>=0.34, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a204=.;
replace a204=_b[_cons] if wave==3;
replace a204=_b[_cons]+_b[onediedtmt] if wave==4;
replace a204=_b[_cons]+_b[onediedtmo] if wave==5;
replace a204=_b[_cons]+_b[onedied] if wave==6;
replace a204=_b[_cons]+_b[onediedtpo] if wave==7;
replace a204=_b[_cons]+_b[onediedtpt] if wave==8;


gen va204=.;

local iter=204;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if evercouple==1 & PI<0.67 & PI>=0.34;

xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI>=0.67, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a205=.;
replace a205=_b[_cons] if wave==3;
replace a205=_b[_cons]+_b[onediedtmt] if wave==4;
replace a205=_b[_cons]+_b[onediedtmo] if wave==5;
replace a205=_b[_cons]+_b[onedied] if wave==6;
replace a205=_b[_cons]+_b[onediedtpo] if wave==7;
replace a205=_b[_cons]+_b[onediedtpt] if wave==8;


gen va205=.;

local iter=205;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets if evercouple==1 & PI>=0.67;

xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI>=0.67, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a206=.;
replace a206=_b[_cons] if wave==3;
replace a206=_b[_cons]+_b[onediedtmt] if wave==4;
replace a206=_b[_cons]+_b[onediedtmo] if wave==5;
replace a206=_b[_cons]+_b[onedied] if wave==6;
replace a206=_b[_cons]+_b[onediedtpo] if wave==7;
replace a206=_b[_cons]+_b[onediedtpt] if wave==8;


gen va206=.;

local iter=206;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if evercouple==1 & PI>=0.67;

*Do medical spending (OOP so use medcost_b) by PI tercile split;

xtreg medcost_b onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.34, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a211=.;
replace a211=_b[_cons] if wave==3;
replace a211=_b[_cons]+_b[onediedtmt] if wave==4;
replace a211=_b[_cons]+_b[onediedtmo] if wave==5;
replace a211=_b[_cons]+_b[onedied] if wave==6;
replace a211=_b[_cons]+_b[onediedtpo] if wave==7;
replace a211=_b[_cons]+_b[onediedtpt] if wave==8;


gen va211=.;

local iter=211;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b if evercouple==1 & PI<0.34;

xtreg medcost_b_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.34, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a212=.;
replace a212=_b[_cons] if wave==3;
replace a212=_b[_cons]+_b[onediedtmt] if wave==4;
replace a212=_b[_cons]+_b[onediedtmo] if wave==5;
replace a212=_b[_cons]+_b[onedied] if wave==6;
replace a212=_b[_cons]+_b[onediedtpo] if wave==7;
replace a212=_b[_cons]+_b[onediedtpt] if wave==8;


gen va212=.;

local iter=212;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b_don if evercouple==1 & PI<0.34;

xtreg medcost_b onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.67 & PI>=0.34, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a213=.;
replace a213=_b[_cons] if wave==3;
replace a213=_b[_cons]+_b[onediedtmt] if wave==4;
replace a213=_b[_cons]+_b[onediedtmo] if wave==5;
replace a213=_b[_cons]+_b[onedied] if wave==6;
replace a213=_b[_cons]+_b[onediedtpo] if wave==7;
replace a213=_b[_cons]+_b[onediedtpt] if wave==8;


gen va213=.;

local iter=213;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b if evercouple==1 & PI<0.67 & PI>=0.34;

xtreg medcost_b_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.67 & PI>=0.34, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a214=.;
replace a214=_b[_cons] if wave==3;
replace a214=_b[_cons]+_b[onediedtmt] if wave==4;
replace a214=_b[_cons]+_b[onediedtmo] if wave==5;
replace a214=_b[_cons]+_b[onedied] if wave==6;
replace a214=_b[_cons]+_b[onediedtpo] if wave==7;
replace a214=_b[_cons]+_b[onediedtpt] if wave==8;


gen va214=.;

local iter=214;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b_don if evercouple==1 & PI<0.67 & PI>=0.34;

xtreg medcost_b onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI>=0.67, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a215=.;
replace a215=_b[_cons] if wave==3;
replace a215=_b[_cons]+_b[onediedtmt] if wave==4;
replace a215=_b[_cons]+_b[onediedtmo] if wave==5;
replace a215=_b[_cons]+_b[onedied] if wave==6;
replace a215=_b[_cons]+_b[onediedtpo] if wave==7;
replace a215=_b[_cons]+_b[onediedtpt] if wave==8;


gen va215=.;

local iter=215;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b if evercouple==1 & PI>=0.67;

xtreg medcost_b_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI>=0.67, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a216=.;
replace a216=_b[_cons] if wave==3;
replace a216=_b[_cons]+_b[onediedtmt] if wave==4;
replace a216=_b[_cons]+_b[onediedtmo] if wave==5;
replace a216=_b[_cons]+_b[onedied] if wave==6;
replace a216=_b[_cons]+_b[onediedtpo] if wave==7;
replace a216=_b[_cons]+_b[onediedtpt] if wave==8;


gen va216=.;

local iter=216;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b_don if evercouple==1 & PI>=0.67;


*Do medical spending (OOP so use medcost_b) by asset split;

xtreg medcost_b onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & below_median_assets_couples==1, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a221=.;
replace a221=_b[_cons] if wave==3;
replace a221=_b[_cons]+_b[onediedtmt] if wave==4;
replace a221=_b[_cons]+_b[onediedtmo] if wave==5;
replace a221=_b[_cons]+_b[onedied] if wave==6;
replace a221=_b[_cons]+_b[onediedtpo] if wave==7;
replace a221=_b[_cons]+_b[onediedtpt] if wave==8;


gen va221=.;

local iter=221;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b if evercouple==1 & below_median_assets_couples==1;

xtreg medcost_b_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & below_median_assets_couples==1, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a222=.;
replace a222=_b[_cons] if wave==3;
replace a222=_b[_cons]+_b[onediedtmt] if wave==4;
replace a222=_b[_cons]+_b[onediedtmo] if wave==5;
replace a222=_b[_cons]+_b[onedied] if wave==6;
replace a222=_b[_cons]+_b[onediedtpo] if wave==7;
replace a222=_b[_cons]+_b[onediedtpt] if wave==8;


gen va222=.;

local iter=222;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b_don if evercouple==1 & below_median_assets_couples==1;

xtreg medcost_b onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & above_median_assets_couples==1, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a223=.;
replace a223=_b[_cons] if wave==3;
replace a223=_b[_cons]+_b[onediedtmt] if wave==4;
replace a223=_b[_cons]+_b[onediedtmo] if wave==5;
replace a223=_b[_cons]+_b[onedied] if wave==6;
replace a223=_b[_cons]+_b[onediedtpo] if wave==7;
replace a223=_b[_cons]+_b[onediedtpt] if wave==8;


gen va223=.;

local iter=223;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b if evercouple==1 & above_median_assets_couples==1;

xtreg medcost_b_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & above_median_assets_couples==1, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a224=.;
replace a224=_b[_cons] if wave==3;
replace a224=_b[_cons]+_b[onediedtmt] if wave==4;
replace a224=_b[_cons]+_b[onediedtmo] if wave==5;
replace a224=_b[_cons]+_b[onedied] if wave==6;
replace a224=_b[_cons]+_b[onediedtpo] if wave==7;
replace a224=_b[_cons]+_b[onediedtpt] if wave==8;


gen va224=.;

local iter=224;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b_don if evercouple==1 & above_median_assets_couples==1;


*Bequests to other heirs:
*gen above_amother=.;
*sum amothers if onediedtmth==1, det;
*replace above_amother=r(mean);
xtreg amothers onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a101=.;
replace a101=_b[_cons] if wave==3;
replace a101=_b[_cons]+_b[onediedtmt] if wave==4;
replace a101=_b[_cons]+_b[onediedtmo] if wave==5;
replace a101=_b[_cons]+_b[onedied] if wave==6;
replace a101=_b[_cons]+_b[onediedtpo] if wave==7;
replace a101=_b[_cons]+_b[onediedtpt] if wave==8;


gen va101=.;
/*replace va101=_se[_cons] if wave==3;
replace va101=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va101=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va101=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va101=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va101=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;

local iter=101;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;


sort wave;
by wave: sum amothers if evercouple==1;

sum amothers if evercouple==1 & onedied==1, d;
sum amothers if amothers>0 & evercouple==1 & onedied==1, d;

*gen above_amother_don=.;
*sum amothers_don if onediedtmth==1, det;
*replace above_amother_don=r(mean);
xtreg amothers_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a102=.;
replace a102=_b[_cons] if wave==3;
replace a102=_b[_cons]+_b[onediedtmt] if wave==4;
replace a102=_b[_cons]+_b[onediedtmo] if wave==5;
replace a102=_b[_cons]+_b[onedied] if wave==6;
replace a102=_b[_cons]+_b[onediedtpo] if wave==7;
replace a102=_b[_cons]+_b[onediedtpt] if wave==8;


gen va102=.;
/*replace va102=_se[_cons] if wave==3;
replace va102=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va102=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va102=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va102=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va102=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=102;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum amothers_don if evercouple==1;


*gen above_amother=.;
*sum amothers if onediedtmth==1 & above_median_assets_couples==1, det;
*replace above_amother=r(mean);
xtreg amothers onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & above_median_assets_couples==1, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a103=.;
replace a103=_b[_cons] if wave==3;
replace a103=_b[_cons]+_b[onediedtmt] if wave==4;
replace a103=_b[_cons]+_b[onediedtmo] if wave==5;
replace a103=_b[_cons]+_b[onedied] if wave==6;
replace a103=_b[_cons]+_b[onediedtpo] if wave==7;
replace a103=_b[_cons]+_b[onediedtpt] if wave==8;


gen va103=.;
/*replace va103=_se[_cons] if wave==3;
replace va103=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va103=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va103=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va103=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va103=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;

local iter=103;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;


sort wave;
by wave: sum amothers if evercouple==1 & above_median_assets_couples==1;

sum amothers if evercouple==1 & onedied==1 & above_median_assets_couples==1, d;
sum amothers if amothers>0 & evercouple==1 & onedied==1 & above_median_assets_couples==1, d;
*gen above_amother_don=.;
*sum amothers_don if onediedtmth==1 & above_median_assets_couples==1, det;
*replace above_amother_don=r(mean);
xtreg amothers_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & above_median_assets_couples==1, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a104=.;
replace a104=_b[_cons] if wave==3;
replace a104=_b[_cons]+_b[onediedtmt] if wave==4;
replace a104=_b[_cons]+_b[onediedtmo] if wave==5;
replace a104=_b[_cons]+_b[onedied] if wave==6;
replace a104=_b[_cons]+_b[onediedtpo] if wave==7;
replace a104=_b[_cons]+_b[onediedtpt] if wave==8;


gen va104=.;
/*replace va104=_se[_cons] if wave==3;
replace va104=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va104=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va104=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va104=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va104=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=104;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum amothers_don if above_median_assets_couples==1;

*gen below_amother=.;
*sum amothers if onediedtmth==1 & below_median_assets_couples==1, det;
*replace below_amother=r(mean);
xtreg amothers onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & below_median_assets_couples==1, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a105=.;
replace a105=_b[_cons] if wave==3;
replace a105=_b[_cons]+_b[onediedtmt] if wave==4;
replace a105=_b[_cons]+_b[onediedtmo] if wave==5;
replace a105=_b[_cons]+_b[onedied] if wave==6;
replace a105=_b[_cons]+_b[onediedtpo] if wave==7;
replace a105=_b[_cons]+_b[onediedtpt] if wave==8;

gen va105=.;
/*replace va105=_se[_cons] if wave==3;
replace va105=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va105=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va105=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va105=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va105=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=105;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

by wave: sum amothers if evercouple==1 & below_median_assets_couples==1;


sum amothers if evercouple==1 & onedied==1 & below_median_assets_couples==1, d;
sum amothers if amothers>0 & evercouple==1 & onedied==1 & below_median_assets_couples==1, d;
*gen below_amother_don=.;
*sum amothers_don if onediedtmth==1 & below_median_assets_couples==1, det;
*replace below_amother_don=r(mean);
xtreg amothers_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & below_median_assets_couples==1, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a106=.;
replace a106=_b[_cons] if wave==3;
replace a106=_b[_cons]+_b[onediedtmt] if wave==4;
replace a106=_b[_cons]+_b[onediedtmo] if wave==5;
replace a106=_b[_cons]+_b[onedied] if wave==6;
replace a106=_b[_cons]+_b[onediedtpo] if wave==7;
replace a106=_b[_cons]+_b[onediedtpt] if wave==8;

gen va106=.;
/*replace va106=_se[_cons] if wave==3;
replace va106=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va106=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va106=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va106=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va106=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=106;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum amothers_don if evercouple==1 & below_median_assets_couples==1;



*gen above_PI_amother=.;
*sum amothers if onediedtmth==1 & evercouple==1 & PI>=0.5, det;
*replace above_PI_amother=r(mean);
xtreg amothers onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI>=0.5, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a107=.;
replace a107=_b[_cons] if wave==3;
replace a107=_b[_cons]+_b[onediedtmt] if wave==4;
replace a107=_b[_cons]+_b[onediedtmo] if wave==5;
replace a107=_b[_cons]+_b[onedied] if wave==6;
replace a107=_b[_cons]+_b[onediedtpo] if wave==7;
replace a107=_b[_cons]+_b[onediedtpt] if wave==8;

gen va107=.;
/*replace va107=_se[_cons] if wave==3;
replace va107=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va107=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va107=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va107=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va107=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/


local iter=107;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum amothers if evercouple==1 & PI>=0.5;

*gen above_PI_amother_don=.;
*sum amothers_don if onediedtmth==1 & evercouple==1 & PI>=0.5, det;
*replace above_PI_amother_don=r(mean);
xtreg amothers_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI>=0.5, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a108=.;
replace a108=_b[_cons] if wave==3;
replace a108=_b[_cons]+_b[onediedtmt] if wave==4;
replace a108=_b[_cons]+_b[onediedtmo] if wave==5;
replace a108=_b[_cons]+_b[onedied] if wave==6;
replace a108=_b[_cons]+_b[onediedtpo] if wave==7;
replace a108=_b[_cons]+_b[onediedtpt] if wave==8;

gen va108=.;
/*replace va108=_se[_cons] if wave==3;
replace va108=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va108=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va108=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va108=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va108=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=108;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum amothers_don if evercouple==1 & PI>=0.5;

*gen below_PI_amother=.;
*sum amothers if onediedtmth==1 & evercouple==1 & PI<0.5, det;
*replace below_PI_amother=r(mean);
xtreg amothers onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.5, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a109=.;
replace a109=_b[_cons] if wave==3;
replace a109=_b[_cons]+_b[onediedtmt] if wave==4;
replace a109=_b[_cons]+_b[onediedtmo] if wave==5;
replace a109=_b[_cons]+_b[onedied] if wave==6;
replace a109=_b[_cons]+_b[onediedtpo] if wave==7;
replace a109=_b[_cons]+_b[onediedtpt] if wave==8;


gen va109=.;
/*replace va109=_se[_cons] if wave==3;
replace va109=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va109=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va109=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va109=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va109=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=109;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum amothers if evercouple==1 & PI<0.5;

*gen below_PI_amother_don=.;
*sum amothers_don if onediedtmth==1 & evercouple==1 & PI<0.5, det;
*replace below_PI_amother_don=r(mean);
xtreg amothers_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & PI<0.5, fe i(HHID);
* generate predicted amother;
* the "constant" is the average fixed effect;
gen a110=.;
replace a110=_b[_cons] if wave==3;
replace a110=_b[_cons]+_b[onediedtmt] if wave==4;
replace a110=_b[_cons]+_b[onediedtmo] if wave==5;
replace a110=_b[_cons]+_b[onedied] if wave==6;
replace a110=_b[_cons]+_b[onediedtpo] if wave==7;
replace a110=_b[_cons]+_b[onediedtpt] if wave==8;


gen va110=.;
/*replace va110=_se[_cons] if wave==3;
replace va110=sqrt(_se[_cons]^2+_se[onediedtmt]^2) if wave==4;
replace va110=sqrt(_se[_cons]^2+_se[onediedtmo]^2) if wave==5;
replace va110=sqrt(_se[_cons]^2+_se[onedied]^2) if wave==6;
replace va110=sqrt(_se[_cons]^2+_se[onediedtpo]^2) if wave==7;
replace va110=sqrt(_se[_cons]^2+_se[onediedtpt]^2) if wave==8;*/;


local iter=110;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum amothers_don if evercouple==1 & PI<0.5;

*Now compute the singles graphs;

local iter = 31;
foreach varname in "assets" "assets_don" "medcost" "medcost_don" "medcost_b" "medcost_b_don" "medicaid_pay" "medicaid_pay_don"{;

sort wave;
by wave: sum `varname';

sort HHID;

*replace prob6=1 if hhstatus==0 & wave==6 & `varname'==.;
* also, make the comparrison variable missing if  
*replace prob6=1 if hhstatus==0 & wave==6 & `iter'==2 & assets==.;
*replace prob6=1 if hhstatus==0 & wave==6 & `iter'==4 & medcost==.;

*egen problem`iter'=sum(prob6), by(HHID);

/*
if use_wave_6==1{;
replace `varname'=. if problem`iter'~=0 & alive5==1;
};
if use_wave_6==0{;
replace `varname'=. if  alive5==1;
};

replace `varname'=. if singlediedtmf==1;
*/
/*
xtreg `varname' singlediedtmo singlediedtmt singledied, fe i(HHID) ;
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a`iter'=.;
replace a`iter'=_b[_cons] if wave==3;
replace a`iter'=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a`iter'=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a`iter'=_b[_cons]+_b[singledied] if wave==6;
*/
*gen b`iter'=.;
*sum `varname' if singlediedtmth==1 & evercouple==0, det;
*replace b`iter'=r(mean);
xtreg `varname' singlediedtmo singlediedtmt singledied if evercouple==0, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a`iter'=.;
replace a`iter'=_b[_cons] if wave==3;
replace a`iter'=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a`iter'=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a`iter'=_b[_cons]+_b[singledied] if wave==6;


gen va`iter'=.;
/*replace va`iter'=_se[_cons] if wave==3;
replace va`iter'=sqrt(_se[_cons]^2+_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[_cons]^2+_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[_cons]^2+_se[singledied]^2) if wave==6;*/;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[singledied]^2) if wave==6;

local iter = `iter' + 1;
};



*gen above_assets=.;
*sum assets if singlediedtmth==1 & above_median_assets_singles==1, det;
*replace above_assets=r(mean);
xtreg assets singlediedtmo singlediedtmt singledied if above_median_assets_singles==1, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a43=.;
replace a43=_b[_cons] if wave==3;
replace a43=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a43=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a43=_b[_cons]+_b[singledied] if wave==6;

gen va43=.;

/*replace va43=_se[_cons] if wave==3;
replace va43=sqrt(_se[_cons]^2+_se[singlediedtmt]^2) if wave==4;
replace va43=sqrt(_se[_cons]^2+_se[singlediedtmo]^2) if wave==5;
replace va43=sqrt(_se[_cons]^2+_se[singledied]^2) if wave==6;*/;
local iter=43;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[singledied]^2) if wave==6;

sort wave;
by wave: sum assets if above_median_assets_singles==1;

*gen above_assets_don=.;
*sum assets_don if singlediedtmth==1 & above_median_assets_singles==1, det;
*replace above_assets_don=r(mean);
xtreg assets_don singlediedtmo singlediedtmt singledied if above_median_assets_singles==1, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a44=.;
replace a44=_b[_cons] if wave==3;
replace a44=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a44=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a44=_b[_cons]+_b[singledied] if wave==6;

gen va44=.;
/*replace va44=_se[_cons] if wave==3;
replace va44=sqrt(_se[_cons]^2+_se[singlediedtmt]^2) if wave==4;
replace va44=sqrt(_se[_cons]^2+_se[singlediedtmo]^2) if wave==5;
replace va44=sqrt(_se[_cons]^2+_se[singledied]^2) if wave==6;*/;
local iter=44;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[singledied]^2) if wave==6;
sort wave;
by wave: sum assets_don if above_median_assets_singles==1;

*gen below_assets=.;
*sum assets if singlediedtmth==1 & below_median_assets_singles==1, det;
*replace below_assets=r(mean);
xtreg assets singlediedtmo singlediedtmt singledied if below_median_assets_singles==1, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a45=.;
replace a45=_b[_cons] if wave==3;
replace a45=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a45=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a45=_b[_cons]+_b[singledied] if wave==6;

gen va45=.;
/*replace va45=_se[_cons] if wave==3;
replace va45=sqrt(_se[_cons]^2+_se[singlediedtmt]^2) if wave==4;
replace va45=sqrt(_se[_cons]^2+_se[singlediedtmo]^2) if wave==5;
replace va45=sqrt(_se[_cons]^2+_se[singledied]^2) if wave==6;*/;

local iter=45;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[singledied]^2) if wave==6;

sort wave;
by wave: sum assets if below_median_assets_singles==1;

*gen below_assets_don=.;
*sum assets_don if singlediedtmth==1 & below_median_assets_singles==1, det;
*replace below_assets_don=r(mean);
xtreg assets_don singlediedtmo singlediedtmt singledied if below_median_assets_singles==1, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a46=.;
replace a46=_b[_cons] if wave==3;
replace a46=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a46=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a46=_b[_cons]+_b[singledied] if wave==6;

gen va46=.;
/*replace va46=_se[_cons] if wave==3;
replace va46=sqrt(_se[_cons]^2+_se[singlediedtmt]^2) if wave==4;
replace va46=sqrt(_se[_cons]^2+_se[singlediedtmo]^2) if wave==5;
replace va46=sqrt(_se[_cons]^2+_se[singledied]^2) if wave==6;*/;

local iter=46;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[singledied]^2) if wave==6;

sort wave;
by wave: sum assets_don if below_median_assets_singles==1;


*gen above_PI_assets=.;
*sum assets if singlediedtmth==1 & evercouple==0 & PI>=0.5, det;
*replace above_PI_assets=r(mean);
xtreg assets singlediedtmo singlediedtmt singledied if evercouple==0 & PI>=0.5, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a47=.;
replace a47=_b[_cons] if wave==3;
replace a47=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a47=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a47=_b[_cons]+_b[singledied] if wave==6;

gen va47=.;
/*replace va47=_se[_cons] if wave==3;
replace va47=sqrt(_se[_cons]^2+_se[singlediedtmt]^2) if wave==4;
replace va47=sqrt(_se[_cons]^2+_se[singlediedtmo]^2) if wave==5;
replace va47=sqrt(_se[_cons]^2+_se[singledied]^2) if wave==6;*/;

local iter=47;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[singledied]^2) if wave==6;

sort wave;
by wave: sum assets if evercouple==0 & PI>=0.5;

*gen above_PI_assets_don=.;
*sum assets_don if singlediedtmth==1 & evercouple==0 & PI>=0.5, det;
*replace above_PI_assets_don=r(mean);
xtreg assets_don singlediedtmo singlediedtmt singledied if evercouple==0 & PI>=0.5, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a48=.;
replace a48=_b[_cons] if wave==3;
replace a48=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a48=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a48=_b[_cons]+_b[singledied] if wave==6;

gen va48=.;
/*replace va48=_se[_cons] if wave==3;
replace va48=sqrt(_se[_cons]^2+_se[singlediedtmt]^2) if wave==4;
replace va48=sqrt(_se[_cons]^2+_se[singlediedtmo]^2) if wave==5;
replace va48=sqrt(_se[_cons]^2+_se[singledied]^2) if wave==6;*/;
local iter=48;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[singledied]^2) if wave==6;

sort wave;
by wave: sum assets_don if evercouple==0 & PI>=0.5;

*gen below_PI_assets=.;
*sum assets if singlediedtmth==1 & evercouple==0 & PI<0.5, det;
*replace below_PI_assets=r(mean);
xtreg assets singlediedtmo singlediedtmt singledied if evercouple==0 & PI<0.5, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a49=.;
replace a49=_b[_cons] if wave==3;
replace a49=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a49=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a49=_b[_cons]+_b[singledied] if wave==6;

gen va49=.;
/*replace va49=_se[_cons] if wave==3;
replace va49=sqrt(_se[_cons]^2+_se[singlediedtmt]^2) if wave==4;
replace va49=sqrt(_se[_cons]^2+_se[singlediedtmo]^2) if wave==5;
replace va49=sqrt(_se[_cons]^2+_se[singledied]^2) if wave==6;*/;

local iter=49;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[singledied]^2) if wave==6;

sort wave;
by wave: sum assets if evercouple==0 & PI<0.5;

*gen below_PI_assets_don=.;
*sum assets_don if singlediedtmth==1 & evercouple==0 & PI<0.5, det;
*replace below_PI_assets_don=r(mean);
xtreg assets_don singlediedtmo singlediedtmt singledied if evercouple==0 & PI<0.5, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a50=.;
replace a50=_b[_cons] if wave==3;
replace a50=_b[_cons]+_b[singlediedtmt] if wave==4;
replace a50=_b[_cons]+_b[singlediedtmo] if wave==5;
replace a50=_b[_cons]+_b[singledied] if wave==6;

gen va50=.;
/*replace va50=_se[_cons] if wave==3;
replace va50=sqrt(_se[_cons]^2+_se[singlediedtmt]^2) if wave==4;
replace va50=sqrt(_se[_cons]^2+_se[singlediedtmo]^2) if wave==5;
replace va50=sqrt(_se[_cons]^2+_se[singledied]^2) if wave==6;*/;
local iter=50;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[singlediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[singlediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[singledied]^2) if wave==6;
sort wave;
by wave: sum assets_don if evercouple==0 & PI<0.5;



*Do this by age  split;

xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age<82, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a401=.;
replace a401=_b[_cons] if wave==3;
replace a401=_b[_cons]+_b[onediedtmt] if wave==4;
replace a401=_b[_cons]+_b[onediedtmo] if wave==5;
replace a401=_b[_cons]+_b[onedied] if wave==6;
replace a401=_b[_cons]+_b[onediedtpo] if wave==7;
replace a401=_b[_cons]+_b[onediedtpt] if wave==8;


gen va401=.;

local iter=401;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets if evercouple==1 & age<82;

xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age<82, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a402=.;
replace a402=_b[_cons] if wave==3;
replace a402=_b[_cons]+_b[onediedtmt] if wave==4;
replace a402=_b[_cons]+_b[onediedtmo] if wave==5;
replace a402=_b[_cons]+_b[onedied] if wave==6;
replace a402=_b[_cons]+_b[onediedtpo] if wave==7;
replace a402=_b[_cons]+_b[onediedtpt] if wave==8;


gen va402=.;

local iter=402;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if evercouple==1 & age<82;

xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age<92 & age>=82, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a403=.;
replace a403=_b[_cons] if wave==3;
replace a403=_b[_cons]+_b[onediedtmt] if wave==4;
replace a403=_b[_cons]+_b[onediedtmo] if wave==5;
replace a403=_b[_cons]+_b[onedied] if wave==6;
replace a403=_b[_cons]+_b[onediedtpo] if wave==7;
replace a403=_b[_cons]+_b[onediedtpt] if wave==8;


gen va403=.;

local iter=403;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets if evercouple==1 & age<92 & age>=82;

xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age<92 & age>=82, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a404=.;
replace a404=_b[_cons] if wave==3;
replace a404=_b[_cons]+_b[onediedtmt] if wave==4;
replace a404=_b[_cons]+_b[onediedtmo] if wave==5;
replace a404=_b[_cons]+_b[onedied] if wave==6;
replace a404=_b[_cons]+_b[onediedtpo] if wave==7;
replace a404=_b[_cons]+_b[onediedtpt] if wave==8;


gen va404=.;

local iter=404;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if evercouple==1 & age<92 & age>=82;

xtreg assets onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age>=92, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a405=.;
replace a405=_b[_cons] if wave==3;
replace a405=_b[_cons]+_b[onediedtmt] if wave==4;
replace a405=_b[_cons]+_b[onediedtmo] if wave==5;
replace a405=_b[_cons]+_b[onedied] if wave==6;
replace a405=_b[_cons]+_b[onediedtpo] if wave==7;
replace a405=_b[_cons]+_b[onediedtpt] if wave==8;


gen va405=.;

local iter=405;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets if evercouple==1 & age>=92;

xtreg assets_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age>=92, fe i(HHID);
* generate predicted assets;
* the "constant" is the average fixed effect;
gen a406=.;
replace a406=_b[_cons] if wave==3;
replace a406=_b[_cons]+_b[onediedtmt] if wave==4;
replace a406=_b[_cons]+_b[onediedtmo] if wave==5;
replace a406=_b[_cons]+_b[onedied] if wave==6;
replace a406=_b[_cons]+_b[onediedtpo] if wave==7;
replace a406=_b[_cons]+_b[onediedtpt] if wave==8;


gen va406=.;

local iter=406;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum assets_don if evercouple==1 & age>=92;

*Do medical spending (OOP so use medcost_b) by age tercile split;

xtreg medcost_b onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age<82, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a411=.;
replace a411=_b[_cons] if wave==3;
replace a411=_b[_cons]+_b[onediedtmt] if wave==4;
replace a411=_b[_cons]+_b[onediedtmo] if wave==5;
replace a411=_b[_cons]+_b[onedied] if wave==6;
replace a411=_b[_cons]+_b[onediedtpo] if wave==7;
replace a411=_b[_cons]+_b[onediedtpt] if wave==8;


gen va411=.;

local iter=411;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b if evercouple==1 & age<82;

xtreg medcost_b_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age<82, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a412=.;
replace a412=_b[_cons] if wave==3;
replace a412=_b[_cons]+_b[onediedtmt] if wave==4;
replace a412=_b[_cons]+_b[onediedtmo] if wave==5;
replace a412=_b[_cons]+_b[onedied] if wave==6;
replace a412=_b[_cons]+_b[onediedtpo] if wave==7;
replace a412=_b[_cons]+_b[onediedtpt] if wave==8;


gen va412=.;

local iter=412;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b_don if evercouple==1 & age<82;

xtreg medcost_b onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age<92 & age>=82, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a413=.;
replace a413=_b[_cons] if wave==3;
replace a413=_b[_cons]+_b[onediedtmt] if wave==4;
replace a413=_b[_cons]+_b[onediedtmo] if wave==5;
replace a413=_b[_cons]+_b[onedied] if wave==6;
replace a413=_b[_cons]+_b[onediedtpo] if wave==7;
replace a413=_b[_cons]+_b[onediedtpt] if wave==8;


gen va413=.;

local iter=413;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b if evercouple==1 & age<92 & age>=82;

xtreg medcost_b_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age<92 & age>=82, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a414=.;
replace a414=_b[_cons] if wave==3;
replace a414=_b[_cons]+_b[onediedtmt] if wave==4;
replace a414=_b[_cons]+_b[onediedtmo] if wave==5;
replace a414=_b[_cons]+_b[onedied] if wave==6;
replace a414=_b[_cons]+_b[onediedtpo] if wave==7;
replace a414=_b[_cons]+_b[onediedtpt] if wave==8;


gen va414=.;

local iter=414;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b_don if evercouple==1 & age<92 & age>=82;

xtreg medcost_b onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age>=92, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a415=.;
replace a415=_b[_cons] if wave==3;
replace a415=_b[_cons]+_b[onediedtmt] if wave==4;
replace a415=_b[_cons]+_b[onediedtmo] if wave==5;
replace a415=_b[_cons]+_b[onedied] if wave==6;
replace a415=_b[_cons]+_b[onediedtpo] if wave==7;
replace a415=_b[_cons]+_b[onediedtpt] if wave==8;


gen va415=.;

local iter=415;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b if evercouple==1 & age>=92;

xtreg medcost_b_don onediedtmo onediedtmt onedied onediedtpo onediedtpt if evercouple==1 & age>=92, fe i(HHID);
* generate predicted medcost_b;
* the "constant" is the average fixed effect;
gen a416=.;
replace a416=_b[_cons] if wave==3;
replace a416=_b[_cons]+_b[onediedtmt] if wave==4;
replace a416=_b[_cons]+_b[onediedtmo] if wave==5;
replace a416=_b[_cons]+_b[onedied] if wave==6;
replace a416=_b[_cons]+_b[onediedtpo] if wave==7;
replace a416=_b[_cons]+_b[onediedtpt] if wave==8;


gen va416=.;

local iter=416;
replace va`iter'=0 if wave==3;
replace va`iter'=sqrt(_se[onediedtmt]^2) if wave==4;
replace va`iter'=sqrt(_se[onediedtmo]^2) if wave==5;
replace va`iter'=sqrt(_se[onedied]^2) if wave==6;
replace va`iter'=sqrt(_se[onediedtpo]^2) if wave==7;
replace va`iter'=sqrt(_se[onediedtpt]^2) if wave==8;

sort wave;
by wave: sum medcost_b_don if evercouple==1 & age>=92;

*1-5 will be couples at death;
*6-10 will be singles at death;
*11-13 will be medical spending;
*21-25 are side bequests;
*30s split by tercile;
*40s split by age group;

gen diff1=-a2+a1;
gen diff2=-a14+a13;
gen diff3=-a16+a15;
gen diff4=-a18+a17;
gen diff5=-a20+a19;

gen sediff1=sqrt(va2^2+va1^2);
gen sediff2=sqrt(va14^2+va13^2);
gen sediff3=sqrt(va16^2+va15^2);
gen sediff4=sqrt(va18^2+va17^2);
gen sediff5=sqrt(va20^2+va19^2);



gen diff6=-a32+a31;
gen diff7=-a44+a43;
gen diff8=-a46+a45;
gen diff9=-a48+a47;
gen diff10=-a50+a49;

gen sediff6=sqrt(va32^2+va31^2);
gen sediff7=sqrt(va44^2+va43^2);
gen sediff8=sqrt(va46^2+va45^2);
gen sediff9=sqrt(va48^2+va47^2);
gen sediff10=sqrt(va50^2+va49^2);

gen diff11=-a4+a3;
gen diff12=-a6+a5;
gen diff13=-a8+a7;

gen sediff11=sqrt(va4^2+va3^2);
gen sediff12=sqrt(va6^2+va5^2);
gen sediff13=sqrt(va7^2+va8^2);


gen diff14=-a34+a33;
gen diff15=-a36+a35;
gen diff16=-a38+a37;

gen sediff14=sqrt(va34^2+va33^2);
gen sediff15=sqrt(va36^2+va35^2);
gen sediff16=sqrt(va37^2+va38^2);


gen diff21=-a102+a101;
gen diff22=-a104+a103;
gen diff23=-a106+a105;
gen diff24=-a108+a107;
gen diff25=-a110+a109;

gen sediff21=sqrt(va102^2+va101^2);
gen sediff22=sqrt(va104^2+va103^2);
gen sediff23=sqrt(va106^2+va105^2);
gen sediff24=sqrt(va108^2+va107^2);
gen sediff25=sqrt(va110^2+va109^2);


gen diff31=-a202+a201;
gen diff32=-a204+a203;
gen diff33=-a206+a205;

gen sediff31=sqrt(va202^2+va201^2);
gen sediff32=sqrt(va204^2+va203^2);
gen sediff33=sqrt(va206^2+va205^2);

gen diff34=-a212+a211;
gen diff35=-a214+a213;
gen diff36=-a216+a215;

gen sediff34=sqrt(va212^2+va211^2);
gen sediff35=sqrt(va214^2+va213^2);
gen sediff36=sqrt(va216^2+va215^2);

gen diff37=-a222+a221;
gen diff38=-a224+a223;

gen sediff37=sqrt(va222^2+va221^2);
gen sediff38=sqrt(va224^2+va223^2);

gen diff41=-a402+a401;
gen diff42=-a404+a403;
gen diff43=-a406+a405;

gen sediff41=sqrt(va402^2+va401^2);
gen sediff42=sqrt(va404^2+va403^2);
gen sediff43=sqrt(va406^2+va405^2);

gen diff44=-a412+a411;
gen diff45=-a414+a413;
gen diff46=-a416+a415;

gen sediff44=sqrt(va412^2+va411^2);
gen sediff45=sqrt(va414^2+va413^2);
gen sediff46=sqrt(va416^2+va415^2);


collapse diff* sediff*, by(wave);
gen SIM=1;

save couples_simulated.dta, replace;
